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ABSTRACT 



Motivated by the limitations encountered with the commonly used direct 
reconstruction techniques of producing mass maps, we have developed a multi- 
resolution maximum-likelihood reconstruction method for producing two dimen- 
sional mass maps using weak gravitational lensing data. To utilize all the shear 
information, we employ an iterative inverse method with a properly selected reg- 
ularization coefficient which fits the deflection potential at the position of each 
galaxy. By producing mass maps with multiple resolutions in the different parts 
of the observed field, we can achieve a comparable level of signal to noise by 
increasing the resolution in regions of higher distortions or regions with an over- 
density of background galaxies. In addition, we are able to better study the 
sub-structure of the massive clusters at a resolution which is not attainable in 
the rest of the observed field. We apply our method to the simulated data and 

OO 

q | to a four square degree field obtained by the Deep Lens Survey. 



Subject headings: cosmology: observations - galaxies: clusters: general - gravita- 
tional lensing 



1. INTRODUCTION 

Recent access to deep, wide field, multi-color optical imaging has established weak grav- 
itational lensing as an effective tool for discovering new clusters of galaxies and measuring 
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the dark matter content of the Universe. Unlike other techniques for selecting clusters, such 
as observing the X-ray emission by the hot intra-cluster medium, the Sunyaev-Zeldovich 
effect on the CMB or over-densities of galaxies in coordinate and/or color space, weak grav- 
itational lensing does not rely on the luminous matter and provides a baryon-independent 
measure of the mass distribution in clusters. 



both on the theoretical front dMiralda-Escudel 



Since the pioneering wo rk of iTvson et al.l (119901). significant progress has been made 



1991 



Kaiser fc Squireslll993l ; ISquires fc Kaiser 



19961 ; ISeitz et al.l Il998l ) and in the quality of the observational data. It is no w possible to 



measure systematic el l ipticities as small as .005 on scales of many arc-minutes (IParker et al. 



20071 ; iJain et al.ll2006l ; IWittman et al.ll2000l ). However, there are still advances to be made in 
reconstructing the surface mass distribution in order to create high resolution and accurate 
mass maps using all of observables such as lens and source redshifts (both spectroscopic and 
photometric), magnification, and distortion measurements from strong lensing arcs. 

In general, there are two classes of reconstruction methods using weak gravitational 
lensing d ata: direct and in v erse m ethods. The direct methods, primarily based on the work 
done by iKaiser fc Squired (119931 ). estimate the surface mass density by approximating a 
local value for the shear (the tidal gravitational field) from the observed ellipticities of the 
background galaxies. In these methods, to avoid the divergence of the statistical uncertainty 
of the surface mass density, the data needs to be smoothed. The smoothing length is a free 
parameter in the reconstruction and there is no a priori way of determining it. 

Furthermore, the transformation of the surface mass density 



K 



an + (1 — a), 



where a is an arbitrary constant does not change the expectation value of the measured ellip- 
ticities ([Schneider &: Seitzlll995l ). Due to this degeneracy in the mass sheet, in the absence of 
redshifts for the sources and lenses, the mass of the clusters cannot be completely determined 
from the shear information alone. In principle, adding extra in formation from other ob serv- 
ables such as magnification can lift the mass sheet degeneracy faroadhurst et al.lll995h . but 
this information cannot be simply incorporated in the direct reconstruction methods. 



Inverse methods aim to find the best fit to the data (ISquires fc Kaiserlll996uBartelmann et al 



19961 ) and have previously only been used to study a small number of individual clusters with 
a low number of background sources. Because of the large number of degrees of freedom, any 
maximum likelihood analysis of weak lensi ng data requires s ome kind o f smoothing to re 



const ruct a mass map. Methods which bin (IBridle et al.lll998l ) or smooth (IBartelmann et al. 



19961 ) the data to regularize the solution do not use all the information provided by each 
source galaxy. To utilize all the shear information, one can employ an iterative inverse 
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method which fits for the deflection potential at the position of each galaxy. Applying the 
smoothing via the a priori expectation of the surface mass density with a properly selected 
regularization weight is a self-consistent, although computationally expensive way of recon- 
structing mass maps from the weak lensing data. Furthermore, the magnification i nformation 



can be simply incorpora ted in inverse methods by constraining the minimization (ISeitz et al. 



1998; Bridle et al. 



1998 



St rong l ensing information from observed arcs can be included in 



the same way (IBradac et al.ll2005l ). 



The necessary uniform smoothing of the data in direct methods of reconstruction (and 
some of the inverse methods) limits the use of the additional information and also produces 
an inconsistent noise level. The resolution of the mass maps is limited by the strength of 
the weak lensing signal and the number density of the background sources which varies 
across the observed field due to Poisson statistics, and also due to background large scale 
structure. In this paper, we present a regularized maximum likelihood method which can 
produce a single mass map with multiple resolutions in the different parts of the observed 
field. Thus, we can achieve a comparable noise level by increasing the resolution at the areas 
with an over-density of sources. Also, the sub-structure of the massive large clusters which 
measurably affect the ellipticities of the background galaxies in a vast area can be studied 
in a map with a resolution which may not be attainable in the rest of the observed field. 
The lower the shear values, the larger the number of sources that are needed to achieve a 
consistent signal to noise per pixel across the filed. Therefore, we can create a map with a 
higher resolution at the areas where the shear values are larger, without reducing the overall 
spatial signal to noise of the detection. 

The only atte mpt to make a multi-resolution mass map is a multi-scale maximum- 
entropy method by Marshall et al.l (120021 ) which uses the intrinsic correlation functions with 
varying width. They report that applying this method to their data did not show a significant 
difference from the single-scale method they also studied. 

The rest of the paper is organized as follows: in § 2, we briefly present the basics of 
gravitational lensing and in § 3, we describe the details of our technique. In § 4, we test the 
method with a simulated one square degree field, distorted with five mock cl usters and in § 5 



we ap ply it to one of the four square degree fields of the Deep Lens Survey (jWittman et al. 
20021 ). Finally, we summarize our method and describe the future plans in § 6. 



-4- 



2. LENSING RELATIONS 

The convergence k, the dimensionless surface mass density of a lens, with a projected 
mass density of S(cc) is defined by fc(cc) = £(:r)/E cr , where S cr = (c 2 /47rG) D s / (D^Dds) 
is the critical surface mass density, and D s are the angular diameter distances of the 
lens and the source from the observer, and is the angular diameter of the lens from the 
source. The angular position in the lens plane is denoted by x. The convergence is related 
to the deflection potential ?p(x), by Poisson's equation k = ^\7 2 ip(x)- The relation between 
the deflection potential and the shear is described by two components, 7 = 71 + 1^2, where 
7i = 3(^,11 - ^,22) and 72 = ^,12, with ip tij = d 2 ij/dx i dx j . 

The images are both distorted and magnified. The distortion is due to the tidal grav- 
itational field, described by shear and the magnification is caused by both isotropic and 
anisotropic focusing, described by both convergence and shear. The ratio of the flux ob- 
served from the image and the true flux of the source defines the magnification, given by 
H= (det^)- 1 = ((1-k) 2 - M 2 )- 1 . 

To study the effects of gravitational lensing on background sources, an ellipse is fitted 
to the shape of each galaxy. The complex ellipticity of a source is defined by the brightness 
quadrupole moments Qij as e = (Qi%— Q22+21Q12) / (Q11+Q22) ■ The transformation between 
the image's quadrupole moments and the source's mo ments due to gravitational 
lensing is given by QW = A Q {s) A kaiser fc Squiresl l 1993), where A is the Jacobi matrix 



A=(1-k) ( l ~ gi ~ 92 V (2) 
V ' \ -92 1 + 9i J 

The parameter g — 7/(1 — k) is the reduced shear. The transformed shapes of objects does 
not depend on the convergence and shear separately, but on the reduced shear. This means 
k and 7 are not direct observables and only their combination in g (or a function of g) can 
be measured from image ellipticities. The shapes are also invariant under the transformation 
g — > 1/g*, so one can define the complex distortion as 



which is invariant under a such transformation (ISchneider fc Seitzlll995l ). 



3. THE METHOD 



In this regularized maximum-likelihood method, our goal is to describe the surface mass 
density distribution over the observed field by fitting for the deflection potential on a multi- 
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resolution grid. We have adapted the prescription of ISeitz et al.l (119981 ) in constructing the 
grid of the deflection potential in order to accommodate for the multi-resolution fitting. In 
a regular maximum-likelihood method, a \ 2 term of the form 



,\2 

X2 N *tZ 



is minimized, where N g is the number of galaxies, is the measured complex ellipticity of 
galaxy k located at (1 < k < N g ), and (e)(a?fc) is the expectation value of the average 
ellipticity at x\. which depends on the deflection potential and its derivatives at that point. 
Should the source redshift information be provided, th e expectation value of th e ellipticity 



could be modified with a cosmological weight function ([Lombard fc Bertinlll999h 



Due to the unmeasurable intrinsic elliptical shape of the weakly distorted galaxies, 
each single object does not provide enough information for lensing reconstruction. Despite 
this, the number of degrees of freedom in minimizing \ 2 is nominally the number of data 
points (i.e. galaxies) subtracted by the number of fitting points. Therefore, the degrees of 
freedom in a wide field reconstruction is quite large and minimizing \ 2 by itself, gives ip 
enough freedom to be able to make it unrealistically small which yields to a potential that 
reconstructs the noise in the data and results in a wrong solution. Adding a regularization 
term with a proper weight helps to constrain \ 2 and avoid over-fitting it. In our method, 
we minimize a function of ip defined as 

T=^ + \n, (5) 

where A is the regularization coefficient. The \ 2 is regularized with a modified zeroth-order 
regularization function 

N X Ny 

^ ^ ^ ^ (^mn Pmn) j (6) 
m=l n=l 

where n mn is the surface mass density at the grid points of the reconstruction grid and p mn 
is the prior. The main advantage of using this function compared to the entropy inspired 
regularization functions is its simplicity. In principle, the entropy inspired functions guaran- 
tee k to be a positive number. But in real numerical analysis when k becomes very small, 
these functions and their derivatives do not behave smoothly. Furthermore, in the presence 
of a mass-sheet degeneracy, it is not clear that enforcing the positivity of kappa is useful in 
deriving a solution. Choosing a function as simple as 1Z which is zero when k is equal to the 
prior is sufficient and ensures the smoothness of the reconstruction. 
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The magnification information can be included by adding to J 7 , a x\ term of the form 



J_ x - (Mi-fi(xi)) 

M 2^ ^ 



v 2 



with A^ number of magnification data points .Mz located at a?; where the predicted value of 
magnification is given by ix(x{). 

The parameter A in equation (jSJ), known as the regularization coefficient, represents the 
compromise between the best fit (i.e the answer that minimizes x 2 ) an d the closest match 
to the prior knowledge (i.e. the answer that minimizes TV). A proper way of finding the best 
value for the regu larization coefficient has been the subject of much debate. For instance, 



Seitz et al.l (119981 ) constrain x 2 to be equal to the degrees of freedom in order to find a good 
guess for A. It is however not clear what a "degree of freedom" is. Because each galaxy shape 
is mostly due to the intrinsic galaxy ellipt i city, t he effective number of degrees of freedom per 



galaxy is much less than 1. iBridle et al.l (119981 ) derive a value for A in a Bayesian manner 



which also has the disadvantage of not easily adapting to the data in a realistic numerical 
analysis. Whatever method employed, the value of the regularization coefficient must be low 
enough so that the solution follows the data and high enough so that it avoids the numerical 
artifacts caused by over-fitting. To determine this value, we minimize T as a function of 
A and compare the resulting x 2 versus 1Z by scaling them to values between and 1. The 
minimum value of x 2 is obtained when A = (the best fit solution, corresponding to Xscaied 
of and TZ sca .i e d of 1), and its maximum value is obtained when only 1Z is minimized (the 
smoothest solution, corresponding to Xscaicd °f 1 an d ^scaled of 0.) The intersection between 
Xscaied vs - ^scaled curve and the line xLaied = "^-scaled determines the proper value of the 
regularizaion coefficient as shown in Figure [H It should be noted that this method dictates 
approximately equal weights to the x 2 term and the regularization term. However, a different 
level of agreement with the prior knowledge (in our case, smoothing) is achieved by selecting 
an intersecting line that has a different slope. Despite the fact that this method requires a 
fair amount of computation time, it ensures the agreement between the data and the a priori 
expectation. 

Our goal is to apply this method to the wide field optical data obtained by the Deep 
Lens Survey. Because the analytical expectations for the average ellipticities of galaxies do 
not take the noise in the data and the shape measurements into account, we estimate the 
expectation value of the average ellipticities as a function of shear based on the simulated 
data. 



In our first suite of simulations, we produce a series of 17.36 square arc-minute simu- 
lated fields in which the simulated galaxies are distributed between magnitudes of 22 and 
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25.5 in the R band. This is approximately the magnitude range of the objects used from the 
DLS data for the mass reconstruction analysis. The ellipticity distribution of the simulated 



galax ies is assumed to be the ellipticity distribution of the galaxies in the UDF (IBeckwith 



20051 ). Because of the small PSF and high signal to noise detection in the UDF data, the 
measured shapes are nearly accurate estimates of the real shapes of the galaxies. Therefore, 
despite the uncertainties due to the finite number of galaxies, the derived ellipticity distri- 
bution is a fair approximation. To include as many galaxies as possible, we assume that 
their average shape does not depend on redshift. We choose the V band data of the UDF to 
determine the ellipticity distribution, because it has the highest signal to noise and is close 
in wavelength to the R band of the Deep Lens Survey, where the shapes of our sources are 
measured. In total, we generate ~ 120,000 simulated galaxies to estimate the expectation 
value of the ellipticity function of distortion. 

We distort the simulated fields according to equation (j5J), varying < g\ < 0.6, while 
#2 is fixed at zero. This distortion step is performed at the pixel-scale of the UDF (0.03 arc- 
seconds). The DLS' PSF is almost always well-sampled, therefore, the simulated images are 
first linearly transformed onto the DLS pixel-scale (0.25 arc-seconds) and then smoothed with 
a Gaussian to simulate the 0.9 arc-second seeing of the data. An appropriate background 
noise is also added to match the simulations to the properties of the actual deep field images. 

Encouragingly, (€2) averages to zero and only (ei) increases with shear in our simulations. 
This demonstrates that it is correct to assume that e\ and €2 have the same orientation as 
gi and g 2 . To avoid the degeneracy between g and l/g*, we use the distortion parameter 
introduced in equation ([3]) and define (e)(xk) = f(\5\ 2 )5. The function /(|5| 2 ) is found from 
the simulations (where g<i = 0, 62 = 0) by 

f(\S\ 2 ) = (e 1 )/5 1 . (8) 
The derived expectation valu es from our simulations along with the analytical approxi- 



mations (ISchneider fc Seitzlll995l ) are shown in Figure [21 Using the same set of simulations, 



we also determine the dispersion in the ellipticity of galaxies a t as a function of distortion. 

Once one exceeds a shear value of 0.6, the shapes of more than 50 per cent of the objects 
in the simulations are not well measured due to splitting by the detection software. This 
causes a bias in the expected ellipticity estimates. Therefore, we extrapolate and use /(|5| 2 ) 
and a e derived from \S\ < 0.88 distortions for the higher values as well. This extrapolation 
is acceptable, because real highly elliptical galaxies are very rare (while high measured ellip- 
ticities are most often caused by unresolved blends of multiple galaxies). Therefore, we can 
and do filter them out of the data without losing any significant weak lensing information. 
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4. IMPLEMENTATION OF THE METHOD 

To create a N x x N y map of the surface mass density distribution, the observed field 
is covered with a (2N X + 4) x (2N y + 4) grid of the deflection potential ipy. Because the 
convergence and shear are second order derivatives of the deflection potential, adding a 
constant or a linear term in x to tp leaves them unchanged and ip needs to be constrained 
to be constant at four of the grid points. For computational simplicity, we have decided to 
keep the four corners of the grid fixed. 

The values of the shear 7^ and convergence at the grid points are obtained by second 
order finite-differencing, hence, the extra rows and columns at each side of the grid. The 
values of shear and convergence at the position of each galaxy in the data are calculated via 
bilinear interpolation. The shear and convergence are computed locally and the coefficients 
relating the deflection potential to n(x k ) and j(x k ) depend only on the geometry of the grid 
and the location of the galaxy at x k , therefore, the coefficients can be calculated once and 
stored to speed up the computations. 

To compute the regularization function 71, we need to know the values of K mn . To fix the 
ringing effects in the projected mass maps caused by second order numerical differentiation 
of ip, we block average ipij at four neighboring grid points and then take the derivatives of 
the deflection potential on this new grid. Hence, the size of the final mass map is N x x N y . 

The components of the shear are also computed on the block averaged deflection po- 
tential grid. Because M.\ are scattered on the field and are not necessarily located on the 
grid points, the matrix H,™ 11 is defined to determine the amount by which each grid point is 
weighted to compute the expected magnification at Xi, changing equation (j7|) to 

1^1/ N *> N * \ 2 

4 = w E- [Mi- E «rvH . (9) 

^ 1=1 M \ m,n=l J 

The elements of 7i only depend on the positions of the magnification data and the structure 
of the grid. Therefore, they too can be calculated once and stored, speeding up the analysis. 
In the presence of magnification data, ip needs to be constrained to be constant only at three 
grid points. 

To min i mize J - ', we use a conjugate-gradient method as encoded in the f rprmn routine by 



Press et al.l (I1992I ). We need to provide this algorithm the first derivatives of T with respect 
to ipij which can be derived with a combination of analytical and numerical methods. In 
general, 

dT dT <97i(^fc) , dT d<y 2 (x k ) dT dK{x k ) 

X — a 1 + a — 7 \ X — a 1 + \ X — q • (10) 



dipij dji(x k ) dijjij d<y 2 (x k ) dK(x k ) dip. 



'.1 
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The derivatives of k, 71 and 72 with respect to ipij only depend on the geometry of the grid 
and the position of the galaxies. They can be derived from the stored coefficients which 
relate the convergence and two components of shear to the deflection potential at the grid 
points. 

The reconstruction procedure starts at a low resolution which depends on the area and 
the number of galaxies of the data. At this level the mass maps are smooth enough and the 
regularization is not required (A = 0). At the end of this step, two very coarse maps of the 
surface mass density and the deflection potential are produced. To increase the resolution, 
we linearly expand and smooth the maps with a Gaussian function (er = 1 pixel, equal to 
the inherent correlation length of the maps) and use the map of ip as the initial potential 
and the map of k as the prior map of the second minimization. 

By finding the proper value of A, we are able to increase the resolution to the limit that 
the data allows us. The resolution of a mass map is limited by the strength of the weak lensing 
signal and the number density of the background sources which varies across the observed 
field due to a variation in source counts and possible background large scale structure. To 
obtain the highest possible resolution, the second step can be repeated: expanding and 
smoothing the tp and k maps of the previous reconstruction and using them as the initial 
potential and prior, respectively. 

In principle, in a maximum likelihood method, the number of unknowns (values of 
the deflection potential on the grid points) must be at least the number of equations (the 
measured ellipticities of the galaxies). However, one single galaxy in the weak lensing limit 
does not provide enough shear information for one grid point because its ellipticity tensor is 
dominated by the random component. Additionally, the Poisson variation in source counts 
and the noise in the shape measurements have undetermined effects on the signal to noise 
across the field. Furthermore, the numerical artifacts in minimizing JF, which also limit 
the resolution are not well predicted. A maximum resolution for a given data set can be 
approximated based on its number of source galaxies, but an exact final resolution of the 
mass maps can not be predetermined. If the signal to noise in a map is not sufficient, we are 
bound to decrease its overall resolution, though we may be able to maintain a high resolution 
at some parts of the field with our multi-resolution reconstruction technique. 

The multi-resolution grid is essentially the same as the single-resolution grid described 
earlier. It only requires an extensive amount of bookkeeping at the edges of the sub-grid 
regions. The shear and convergence computations for the galaxies in the middle regions of 
the sub-grids are performed similarly to the single-resolution computations. For the galaxies 
which lie on the edge or corner cells, the values of the deflection potential at the required 
positions in the field with no real grid points allocated for them are interpolated. As in 
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the single-resolution construction, the coefficients relating the convergence and shear to the 
deflection potential depend only on the position of the galaxies and the geometry of the main 
grid and the sub-grids, thus this step is required to be performed only once. 

In order to simplify the calculations, the resolutions of the rectangular sub-grids, which 
may be different from one to another, are required to be 2n times higher than the original 
resolution. The maps of ij) and k produced in the final single-resolution reconstruction are 
used as the initial potential and the prior, respectively. The proper regularization coefficient 
is derived similarly to the single-resolution reconstruction. The minimization of the function 
T is performed over all grid points in the main grid and sub-grids, except for the four corners 
that are held constant. 



5. SIMULATED DATA 

We simulate a one square degree field distorted by 5 clusters with the NFW profile 
(jNavarro. Frenk fc White! 119971 ) at a redshift of 0.4 with masses ranging between 10 13 to 



10 15 Solar masses. The mass and position of each cluster is detailed in Table [T] and the 
analytical expectations of the surface mass density map due to these clusters is shown in 
Figure [3] (top left). Clusters number 1 and 2 are chosen to be close to each other to test our 
ability to separate bright adjacent peaks using the multi-resolution method. Cluster number 
3 is a typical isolated cluster and clusters 4 and 5 are are intentionally chosen to be low-mass 
clusters to study the lower signal to noise limits of the reconstruction by our technique. 

The angular diameter distances are evaluated assuming a ACDM cosmology with Q m = 
0.3, — 0.7, and the Hubble constant H = 70 km s _1 Mpc -1 . The objects are randomly 
oriented galaxies. The ellipticity distribution is assumed to be the ellipticity distribu tion of 



the g alaxies in the UDF and their number density follows a power law distribution (ITyson 



19881 ). The galaxies are divided between seven redshift layers based on their magnitudes, 
which range between 23 and 27 in the R band: Z\ = 0.27, z 2 = 0.45, z 3 = 0.68, z± = 0.90, 
z 5 = 1.20, zq = 1.95 and z-j = 3.00. These logarithmically determined layers fairly simulate 
the redshift distribution of the galaxies in the Deep Lens Survey and varying these values, 
especially the furthest redshift, does not change the total distortion by a measurable amount. 
After distorting, we convolve the image to a seeing of 0.9 arc-seconds. 

To measure the sh ape of the galaxies, w e employ the same procedur e used in the Deep 



Lens Survey's pipeline (jWittman et al.ll2006l ). Briefly, we use SExtractor (IBertin fc Arnouts 



19961 ) to detect the objects. The improve upon the shape measurements which are not 



optimal for weak lensing studies, we employ the ellipto program, which can produce more 
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accurate shape measurem ents via an iterative weigh ting algorithm, where the weight function 



is an elliptical Gaussian (IBernstein &: Jarvisll2002l ). We apply the same selection criteria in 



magnit ude and size applied to the DLS data to select objects to be used in making the mass 



maps (jWittman et al.ll2006l ). We require that t he moments be su c cessfu lly measured by 



ellipto and employ the size measure defined by IBernstein fc Jarvid (120021 ) to filter out the 
objects smaller than the PSF (ellipto-size of 5 pixel 2 ). We only keep the objects brighter 
than the magnitude 25.5. After filtering out the unwanted objects, we have a catalog of 
~ 109,000 galaxies. 

We start off the reconstruction at a resolution of 3 arc-minutes per pixel on the grid of 
the deflection potential with a constant initial value over the field which yields a mass map 
with a 6 arc-minute per pixel resolution. At this level, the maps are coarse enough that there 
is no need for any regularization, hence A = 0. To find the proper regularization coefficient 
for the higher resolution reconstructions we follow our recipe and run minimizations with 
coefficients between A = and A = 10 in addition to minimizing only 1Z at each step to 
finally produce a 1'.5 per pixel mass map. Figures [1] and [3] (top right) show xLaied vs - ^scaled 
for the last step of this reconstruction process and the final mass map, respectively. Although 
we do not probe the entire parameter space directly at the highest resolution, we vary the 
values of the deflection potential evenly over the lowest resolution grid with small and large 
increments which does not produce a lower x 2 , assuring that the conjugate-gradient method 
reaches the minimum and does not stop at a possible local minima. 

Due to the low signal to noise detection of the lowest mass clusters, it is not possible 
to increase the resolution of the overall map. However, it is still possible to increase the 
resolution at the vicinity of the first and second clusters, where we increase the resolution of 
the mass map in a square region by a factor of two to 0'.75 per pixel. In the single-resolution 
map, these clusters are reconstructed without any separation (i.e. as a single object). The 
resulting multi-resolution convergence map (Fig. [31 bottom right) shows the cluster not only 
with the expected symmetric profile, but also very well separated (with the peaks detected 
at 2'. 65 of each other, in very good agreement with the 2'. 9 separation of the input profile). 

Because of the differential nature of our fitting function, the pixels of the mass maps 
created by our method are not strongly correlated with each other. Therefore, the total 
surface mass density of each deflector can be measured by summing over the values of the 
pixels which are above a predetermined threshold. We measure Ktotai within the r2oo radius 
of each deflector, setting the detection threshold at 2 times the background rms. At this 
threshold level, all five deflectors along with three spurious objects are detected. Obviously, 
increasing the detection threshold will remove the spurious objects, however, the weakest 
deflector would not be detected either (for instance at 3 times the background rms). In 
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the absence of other observational data such as redshift or magnification information, the 
mass sheet degeneracy cannot be broken. Nonetheless, our measurements (Table [2]) are in 
close agreement in positions and total surface mass densities w ith the measurements from an 
analytically calculated map of convergence (Fig. [31 top left) (jWright fc Brainerdl 120001 ) . A 
mass sheet corresponding to the degeneracy coefficient of a ~ 0.88 (Eqn. [1]) transforms the 
measured surface mass density to the expected surface mass within the estimated errors. The 
effects of this degeneracy in our inverse method are most probably suppressed, because the 
reconstruction process is started with the assumption that the field is empty of any structure. 
This is an initial condition that cannot be incorporated in a direct method reconstruction. 

We also reconstruct the convergence map of a catalog made by distorting the same 
simulated source galaxies with five deflectors located at the same position but with half the 
stren gth (i.e. the M 2 oo mass of each cluster is reduced by half.) As expected (IBridle et al. 



19981 ). the noise level and the regularization coefficient at each step of the reconstruction 
remain the same as the original reconstruction process. However, the two weakest deflectors 
are not detected at all when the detection threshold is set at 3 times the background rms. We 
similarly reconstruct the convergence map of a catalog distorted by the original deflectors but 
with only half the background galaxies. The change in the number of sources also changes 
the regularization coefficient. After determining the proper value of A and making the final 
mass map, the measured signal of the three more massive clusters is very close to the signal 
measured from the original mass map while the two least massive ones are not detected. 



I n addition, we reconst ruct the surface mass density employing a direct method (IKaiser fc Squires 



19931 ; IWittman et al.ll2006l ). using the weight function introduced by lFischer fc Tysonl (119971 ) 



W(r) = (1 - e^( r2 Af n )) e -|(r 2 /^ut) j 



(11) 



with r in = l'.l and r out = 12'. 5. The atmospheric and optical distortions of the shapes of the 
background sources res ult in suppressed sig n als. W e correct for these effects by employing the 
method introduced by iBernstein fc Jarvid (120021 ) and approximate the amount of required 
adjustments to the ellipticities of each source galaxy. In the resulting mass map (Fig. [31 
bottom left), when the detection threshold is set at 2 times the background rms, we are able 
to detect all five deflectors along with nine spurious objects. 

The pixels in the direct method map are highly correlated. Moreover, because of the 
weight function (Eqn. [TT]) . it is the convolved surface mass density that is measured from 
this map. Therefore, it is not proper to compare the Ktotai measurements with the previous 
measurements, and thus the direct reconstruction map is only suitable to study the number 
count of clusters and possibly the relative strength of their signal. 



To estimate the statistical significance of detecting clusters at different resolutions given 
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our data, we perform a set of Monte Carlo simulations and create a number of source catalogs 
in which the ellipticity components of one galaxy is given to another, though their positions 
are not changed. The mass map for each catalog is created by starting at the initial resolution 
of the original mass map and the same procedure is followed to achieve the final resolution 
using the same regularization coefficients of the original reconstruction process at each step. 
Figure H] shows that there are not any objects in the Monte Carlo catalogs with signals larger 
than or equal to the combined signal of the first and second clusters, where the detection 
threshold is set at 1.5 times the background. This is also true for the third cluster. The 
histogram in Figure H] can also be interpreted as the probability distribution that the peaks 
are real detections. We calculate the probability of measuring a signal within the r 20 o radius 
of each deflector that is equal to its fttotai by measuring the probability of finding the same 
signal in randomly selected regions of the Monte Carlo mass maps (Table [2]). When there are 
no detected objects with a given signal, a rough lower limit for the probability of detection 
being real can be estimated by the inverse of the number of th e Monte Carlo simula tions 
per detected objects in the original catalog with that signal (IWall fc Jenkins! 120031 . and 
references within). In addition, because the r 2 oo of each cluster is a known priori, we can 
estimate the l-a error for the measured total surface mass density of the clusters, using the 
same set of Monte Carlo simulations. This Monte Carlo analysis shows that we have been 
able to detect the more massive clusters with a high probability of being real detections and 
also measure their total surface mass density in good agreement with the analytical input. 
The total surface mass density measurements for the lower mass clusters are also in good 
agreement with the analytical input. However, the high number of detected objects in the 
Monte Carlo simulation with similar signals to those of the less massive clusters, suggests a 
lower probability that any detection peak is a real object. 



6. WIDE FIELD OPTICAL DATA 

We also apply our method to reconstruct the mass distribution over a 4 deg 2 field with 
deep optical imaging (R < 26), obtained by the Deep Lens Survey. The DLS is a multi-color 
survey of five separate patches of sky with a consistently good image quality (< 0.9") in the 
R band (where the shapes of the source galaxies are measured). We do not intend to break 
the mass sheet degeneracy in this paper and only use the shear information in the data. We 
run our method on the DLS field 2 (F2) centered at RA = 09 h 19 m 32 s .4, DE C = +30 o 00'00". 



For t he weak lensing analysis, the data is cleaned of unsuitable objects (IWittman et al. 



20061 ). Stars and any object smaller than the PSF size are removed, using the ellipto-size 
vs. magnitude diagram. The bright end of the locus which contains saturated objects and 
bright galaxies is also filtered out. We also only keep the galaxies with successfully measured 
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intensity moments (by ellipto) which are brighter than R = 25 to reduce the noise due 
to the faintest and noisiest galaxies. After filtering the unwanted objects out, there are 
~ 140,000 galaxies left in the data set (Fig. ED- 

In the same way as described in the previous section, we start the reconstruction process 
at a very low resolution of 6 arc-minutes per pixel without regularizing the x 2 > that produces 
a 12 arc-minute per pixel mass map. The process is continued and the higher resolution mass 
maps with the appropriate regularization coefficients are created. After four steps, the final 
mass map with a resolution of 1'.5 per pixel is created (Fig. El left). This figure (right) also 
shows the direct mass reconstruction of this field with r- m = 2'. 9 and r out = 24'. 4 (Eqn. [TT]) . 

The largest signal in this field is due to a set of known clusters (the Abell 7 81 complex) 



which consists of several independent components at redshifts of 0.29-0.43 (IGeller et al. 



20051 ). In the final single- resolution mass map, the sub-structure of this system is not very 
well resolved. However, the signal due to this complex is high enough to allow a higher 
resolution reconstruction which the rest of the field does not permit. Therefore, an area 
(0.09 deg 2 ) around this region for the multi-resolution reconstruction is chosen. The resulting 
mass map is shown in Figures [7] and El in which three out of the four spectroscopically 
confirmed components of this system are very well resolved. Two other bright peaks also 
appear in the vicinity of this system, which will require more investigation to be confirmed. 
We also perform the multi-resolution reconstruction on two random regions of this field void 
of areas with large signal. The result is mass maps in which the noise has been fitted for 
rather than the signal, showing that a higher global resolution is not attainable with this 
source catalog (Fig. [9]). 

The same Monte Carlo method described earlier is employed to estimate the statistical 
significance of detecting clusters in this field. Neither the r 2 oo radii nor the redshifts of the 
cluster candidates in this field are a priori known. Therefore, we measure the total isophotal 
signal, setting the detection threshold is set at 1.5 times the background rms. Figure [TOl 
shows the number of the detected objects with a given total signal per catalog. This graph 
indicates that the number of detected objects per catalog with signals larger than or equal 
to those of the top two cluster candidates is insignificant, thus they are detected with very 
high signal to noise and their realness is highly probable. However, the high number of 
objects per one Monte Carlo catalog with signals equal to the lower ranking objects in the 
DLS field suggests that these objects have a much lower probability of being real detections. 
Conversely, the results implifies that a significant number of "clusters" detected at this level 
are spurious. 
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7. CONCLUSION 

In this paper we have introduced a maximum-likelihood method for weak lensing con- 
vergence map rec onstructions. This method, which is primarily based on the prescription of 



Seitz et al.l (119981 ) is able to produce multi-resolution mass maps that can be used to achieve 
comparable noise levels in regions of higher distortion or regions with an over-density of 
background sources. In addition, the sub-structure of massive clusters can be better studied 
at a resolution that is not attainable in the rest of the field. The expectation value of the 
ellipticities of sources is estimated via realistic simulations and the regularization coefficient 
is properly chosen to be what the data dictates itself. 

We test the performance of our method on a one square degree simulated field and 
conclude that reconstructing mass maps does not depend on the initial conditions. Although 
we did not expect to break the mass sheet degeneracy, our surface mass density measurements 
are in good agreement with the analytical expectation. The effects of this degeneracy seem 
to be suppressed in the simulations, because the reconstruction process is initiated with the a 
priori assumption that there are no structures in the field. The relatively high source number 
density of the simulated field (~ 30 galaxies per square arc- minute), is only sufficient to 
detect the top four massive deflectors with high signal to noise and the fifth ranking cluster 
(M 2 oo ~ 0.7 x 10 14 Solar masses) is not detected when the detection threshold is set to 
remove all spurious detections. Reducing the source number density to ~ 15 galaxies per 
square arc-minute, lowers the signal to noise for the less massive clusters and both fourth 
(M200 ~ 1-3 x 10 14 Solar masses) and fifth ranking clusters are not resolved. However, the 
total surface mass density of the top three clusters measured from the low source density 
catalog is very similar to the previous measurements from the original catalog. In addition, 
we reconstruct a multi-resolution mass map of this field with the highest resolution of 0'.75 
per pixel, in which the first and second clusters are successfully separated and the expected 
symmetric profiles are resolved. The Monte Carlo type simulations created by shuffling the 
ellipticities of the source galaxies in the simulated field demonstrate that the less massive the 
clusters, the higher the number of detected objects with similar signal, solely due to random 
orientation of background sources. From these simulation, we also estimate the probability 
for the peaks' detections to be real. 

We also report a preliminary convergence map of a 4deg 2 field obtained by the DLS 
and reconstruct a multi-resolution mass map. This map, unlike the single-resolution one, 
successfully shows the sub-structure of the brightest system in the field, corresponding to 
the Abell 781 complex, clearly resolving three of its components. Employing Monte Carlo 
simulations, we show that only the top two cluster candidates in the single-resolution map 
have a significant probability of being real clusters whereas the realness of the rest of the 
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candidates is not highly probable. 

Mass reconstruction by this multi-resolution inverse method can be improved in many 
ways. The redshift information of the background sources can be easily incorporated in 
the expected ellipticity function. This method is also capable of including other available 
observational information such as magnification data in the lensing reconstruction. The 
application of this method to the DLS data set will be the first attempt in breaking the 
degeneracy in wide field mass reconstruction using both shear and magnification data. Papers 
presenting the mass function and the biases in the mass reconstruction of this field with a 
more comprehensive analysis of the confirmed shear selected clusters, as well as the statistical 
properties of candidate systems are in preparation. 
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Fig. 1. — The resulting xLaicd vs - "^scaled from 9 minimizations with coefficients between 
A = and A = 10 in addition to the result from minimizing only TZ. The intersection 
between the Xscaicd^scaicd curve and the line xLaied = ^scaled indicates A = 0.001 to be the 
proper coefficient for the 1'.5 per pixel reconstruction of the simulated field. 
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Fig. 2. — The expectation value of e\ and e 2 versus \8\ 2 . Solid squares represent the simula- 
tions at the DLS pixel-scale with simulated PSF and matched signal to noise, open squares 
represent the simulations of a UDF like field and s olid triangles represent the analytically 
approximated expectations ( Schneider fc Seit3 1995 ). 



Fig. 3. — The input surface mass density of the five clusters listed in Tabled] (top left) and 
the single-resolution inverse reconstruction (top right), the direct reconstruction (bottom 
left) and the multi-resolution inverse reconstruction (bottom right) of the 1 deg 2 simulated 
field. The resolution in the box around the first and second clusters in the multi-resolution 
map is 0'.75 per pixel whereas the resolution in the rest of the field is 1/.5 per pixel (~ 480 
Kpc at z — 0.4), which is also the resolution of the other maps. 
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Fig. 4. — The histogram of the number of detected objects ti(k) per catalog for a given total 
surface mass density «totai- The solid histogram shows n(/c) for the Monte Carlo catalogs 
base on the simulated field catalog (made by shuffling the ellipticities of the source galaxies in 
the simulated field) and the dashed histogram shows n(n) for the original simulated catalog. 




Fig. 5. — The ellipto-size vs. magnitude diagram of the objects in a random 0.1 deg 
region of the DLS F2. To clean the data of unsuitable objects for the weak lensing analysis, 
stars and any object smaller than the PSF size, along with the bright end of the locus 
which contains saturated objects and bright galaxies are filtered out. We also only keep the 
galaxies with successfully measured intensity moments (by ellipto) that are brighter than 
R = 25, selecting o nly th e objects inside the box to be consistent with the selection method 
of Iwittman et all tood ). 




Fig. 6. — The final single-resolution convergence maps of the 4 deg 2 DLS F2, from the inverse 
(left) and the direct (right) reconstructions, with a resolution of 1'.5 per pixel. 
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Fig. 7. — The final multi-resolution convergence map of the 4deg 2 DLS F2. The resolution 
in the box around the A781 complex is 0'.75 per pixel whereas the resolution in the rest of 
the field is 1'.5 per pixel. 
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Fig. 8. — The section of the final multi-resolution convergence map of the DLS F2 (Fig. [7]) 
in the vicinity of the the A781 complex with overlayed intensity contours. The resolution in 
the box around the A781 complex is 0'.75 per pixel whereas the resolution in the rest of the 
field is V.5 per pixel. 
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Fig. 9. — The sections of the multi-resolution reconstruction of a random region (left) and 
the original single-resolution reconstruction of the same region (right) of the DLS F2. This 
region is void of areas with large signal. The resolution in the box is 0'.75 per pixel whereas 
the resolution in the rest of the field is 1'.5 per pixel. It is clear that in the multi-resolution 
region, it is the noise that has been fitted for rather than the signal, indicating that a higher 
global resolution is not attainable with this source catalog. 
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Fig. 10. — The histogram of the number of detected objects n(n) per catalog for a given total 
surface mass density «totai- The solid histogram shows n(/t) for the Monte Carlo catalogs 
base on the DLS field catalog (made by shuffling the ellipticities of the source galaxies in the 
simulated field) and the dashed histogram shows ti(k) for the original DLS F2 catalog. 
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Table I. Simulated NFW Clusters 



Cluster 


x (pix) 


y (pix) 


r a (Kpc) 


T-200 (Mpc) 


Mass (10 14 x M Q ) 


1 


10000.0 


10000.0 


430.13 


2.151 


26.1 


2 


9500.0 


9500. 


268.83 


1.075 


3.3 


3 


4000.0 


3500.0 


322.60 


1.505 


9.0 


4 


5000.0 


8000.0 


172.05 


0.806 


1.3 


5 


8250.0 


5400.0 


134.42 


0.645 


0.7 



Note. — Properties of the simulated NFW clusters (z = 0.4). The height and 
width of the field are 1 degree = 14400 pixels. 
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Table 2. Measured fttotai of Simulated Clusters 



Analytical Input Inverse Method 



Cluster 


Position (pix) 


K total 


Position (pix) 


K total 


-frcal 


1, 2 


(26.76, 26.76) 


5.310 


(26.87, 27.43) 


6.164 ± 0.360 


99.97% 


3 


(11.03, 9.70) 


1.510 


(11.15, 9.66) 


1.640 ± 0.273 


99.95% 


4 


(13.76, 21.80) 


0.308 


(14.00, 22.44) 


0.250 ± 0.129 


84.83% 


5 


(22.22, 14.79) 


0.170 


(23.52, 15.00) 


0.134 ± 0.096 


72.62% 



Note. — The measured total surface mass density of the simulated clusters from 
the analytical input and our inverse method, all shown in Figure [3] A mass sheet 
corresponding to the degeneracy coefficient of a ~ 0.88 (Eqn. [TJ transforms the 
measured surface mass density to the expected surface mass within the estimated 
errors. The error and probability estimates are derived from Monte Carlo simulations. 
The probability of finding objects in randomly selected regions of the Monte Carlo 
mass maps with the same or less signal than that of each cluster determines the 
probability of detecting such signal solely due to random orientation of background 
sources. One minus this probability is a fair estimate for the probability of detections 
to be real, P rea i- 



